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ABSTRACT 

The key to using a strong gravitational lens system to measure the Hubble constant 
is to obtain an accurate model of the lens potential. In this paper, we investigate the 
properties of gravitational lens B1608+656, a quadruply-imaged lens system with an 
extended source intensity distribution. Our analysis is valid for generic quadruply- 
lensed systems. Limit curves and isophotal separatrices are defined for such systems, 
and we show that the isophotal separatrices must intersect at the critical curves and 
the satellite isophotes must be ta ngent to the limit curves. The current model of 
B1608+656 l|Koopmans et a l. 2003) satisfies these criteria for some, but not all, of the 
isophotal separatrices within the observational unc ertainty. We study a non-parametric 
method of potential reconstruction proposed by iBlandford. Surpi fc Kundicl (2001) 
and demonstrate that although the method works in principle and elucidates image 
formation, the initial potential only converges to the true model when it is within ~ 1 
percent of the true model. 
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1 INTRODUCTION 

Strong gravitational lens systems provide a tool for mea- 
suring cosmological parameters. With the measured relative 
arrival time delays between the multiple images of the lensed 
source and a model of the lens potential, one can deduce a 
value of the Hubble constant. In addition, strong gravita- 
tional lens systems can be used to probe galaxy mass dis- 
tributions, including dark matter, since the le ns potential is 
directly related to the lens mass distribution. llRefsdall964l) 
Several strong gravitational lenses with either two im- 
ages ("doubles") or four images ("quads") have been ob- 
served. The ones with extended source distributions are of 
special interest since they provide additional constraints for 
the lens potential due to surface brightness conservation. 
The traditional approach to modelling the lens mass dis- 
tribution is to postulate a parametric form for the lens 
distribution and minimize some chi-square to fit the data. 
The method is limited by the choice of the parameters; 
as the observational quality improves, the original para- 
metric model generally b ecomes inadequate to fit the data 
JWilliams fc Sahal (2000) consider a pixellated mass distri- 
bution which is non-parametric, but use only the nuclear im- 
age positions and not information from the extended source 
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to constrain the distribution). Ideally, we want a method 
that employs the extended source information to obtain a 
non-parametric form of the lens potential whose accuracy 
is limited solely by the observational noise in the data. 
iKoopmansI i2005l) has also taken this approach. 

In section 2, we study in detail a quadruply im- 
aged gravitational lens system with an extended source, 
B1608+656, showing the various criteria that the isophotes 
of the extended source must satisfy. In section 3, we ex- 
amine a method of p otential reconstruction proposed by 
IBlandford et alJ i200ll) to correct the potential values pixel 
by pixel from a starting perturbed potential model. 



2 PROPERTIES OF A QUADRUPLY LENSED 
SYSTEM 

We will focus on the quadruply imaged gravitational lens 
system B1608+656 in this section. Fig. Q shows an image 
of the system take n by the Hubble Space Te lescope through 
the F814W filter JSurpi fc Blandfordll2003l). The source is 
at a redshift of z s = 1.39 iFassnacht et alJll996l) and its 
images are labelled by A, B, C, and D. The system has two 
lens galaxies, Gl (the primary lens) and G2 (the secondary 
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Figure 1. The original reduced ffi?T/F814W image of 
B1608+656. The four images are labelled A, B, C, and D ; the 
two lens galaxies arc Gl and G2. iSurpi fc Blandfordll2003t simi- 
lar image in lKoopmans et alj|200a) 



lens), that are at a redshift of Zd = 0.63 iMvers et aliWddh . 1 
The galaxy Gl is about five times more massive than G2. 
B1608+656 is unique in that all three relative time delays 
between the four images are determined with accuracies of a 
few percent. The time delays relative to image B for images 
A, G and D are 31. 5±1.5 days, 36.0±1.5 days and 77.0±1.5 
days, respectively dFassnacht et aljll999l . l20o3) . 

Section 2.1 that follows is a review of the theory of 
gravitational lensing. Readers familiar with lensing may 
wish to proceed directl y to Section 2.2, which analyses the 
iKoonmans et all l)2003h model of B1608+656. 



In terms of the dimensionless surface mass density, denoted 
by k(6), the lens potential is 



- [ d 2 9'k(8') in \e-e'\ 



(3) 



The time delay function relative to the case of no lensing 
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where Dd, D s , and Dd s are, respectively, the angular diam- 
eter distance from us to the lens, from us to the source, and 
from the lens to the source. 

The constant coefficient in equation JIJ is proportional 
to the angular diameter distance and hence inversely pro- 
portional to the Hubble constant in a flat A-CDM universe. 
Therefore, by measuring the relative time delays between 
the various images, we can in principle deduce the value of 
the Hubble constant if we know the source position ((3) and 
the lens potential (ip(0)). 

To characterise the magnifications of images in gravita- 
tional lensing, a Hessian is used 



(5) 



Using the lens equation Q , the above equation can be writ- 
ten as 



A(0) 
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where the subscript 1 (or 2) in ip indicates a derivative with 
respect to 0\ (or 82). The magnification matrix is defined as 
fx = A -1 , and the associated magnification factor is 



(7) 



det A(0) 

According to equation Q, the positions with det A(0) = 
have divergent magnification; the loci of such points on the 
image plane define the critical curves. Using the lens equa- 
tion critical curves on the image plane are mapped to 
caustic curves (or simply caustics) on the source plane. The 
caustic curves separate regions of different image multiplic- 
ities. 



2.1 Gravitational lensing 

Readers familiar with gravitational lensing may wish to skip 
this se ction. We follow lKochanek. Schneider fc Wambsgansa 
i2004l) for the theory of gravitational lensing. 

Let us denote the angular coordinates on the source and 
image planes by (3 — (Pi,^) and = (81,62), respectively. 
The lens equation governing the deflection of light rays is 



(3 = - a{6), 



(1) 



where a(G) is the scaled deflection angle that is the gradient 
of a scalar function called the lens (or deflection) potential: 



a(e) = vy>(0). 



(2) 



The quoted redshift is that of Gl. We assume that G2 is at the 
same redshift as Gl since G2 is too faint for its redshift to be 
measured. 



2.2 Gravitational lens B1608+656 

To investigate the anatomy of the quad B 1608+656, we use 
the m ass distribution model proposed by iKoopmans et all 
(2003). The parametric form of the dimensionless surface 
mass density for each of the two lens galaxies is a singular- 
power law ellipsoid: 
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where (0 g ai l ,dgai 2 ) are coordinates relative to the galaxy 
centre and b, qi, and 7' are parameters to fit the data. The 
origin of coordinate 6 is set at the position of image A. Each 
of the lens galaxies is centred at the coordinates (On, 612) and 
is rotated by a major-axis position angle Op a that is mea- 
sured from north to east (top to left). There is an additional 
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Table 1. Values for parameters in equations |HJ and JHJ for 

the B 1608+656 SPLEl+D(isotropic) model in Koopma ns et alj 
<200fl 



lens galaxy Gl G2 

b 0.526 0.269 

qi 0.604 0.318 

■y' 2.05 2.12 
centroid (0n,0I2) (0.425,-1.069) (-0.291,-0.928) 

position angle 6 PA (°) 77.0 68.4 

"text 0.077 

shear position angle 9 ex t{°) 13.4 



external shear centred on Gl whose contribution to the lens- 
ing potential, in polar coordinates relative to the shear cen- 
tre ((r,cj>) with 8 s h 1 = rcos((f>) and s h 2 = rsm.((j>)), is 

1 2 

Ipext(Osh) = ^extT COS(20), (9) 

where ^^xt is a parameter characterising the shear strength. 
The rotation of the external shear is given by the po- 
sition angle 8 e xt- We adopt the parameter va l ues o f the 
SPLEl+D(isotropic) model in lKoopmans et alJ i200ot) and 
list them in Tabled 



2.2.1 Critical and caustic curves 

The critical curves on the image plane and the caustic curves 
on the source plane of the SPLE1+D (isotropic) model in 
iKoopmans et all J2003f) are shown in Fig. [5] in the middle 
panel and the left panel, respectively. The locations of the 
lens galaxies are indicated by open triangles on the image 
plane. The marked source and image locations will be dis- 
cussed in the next section. With the two elliptical lens galax- 
ies, the large critical curve loop is a deformed version of an 
elliptical curve of one singular power law ellipsoid (equa- 
tion JSJl). The corresponding diamond shaped caustic curve, 
known as an astroid, is typical for elliptical mass distri- 
butions. An astroid is composed of four folds (branches of 
smooth curves) joining at four cusps. An individual power 
law ellipsoid has an astroid that is symmetrical with respect 
to the semi-major and semi-minor axis of the lens. With 
the two lens galaxies in the SPLEl+D(isotropic) model, we 
have an asymmetry in the astroid and two additional small 
triangular caustics, called the deltoids, that map into the 
small loops on the image plane. 



2.2.2 Image positions and time delay surface 

It is instructive to see how the images move on the image 
plane as the source is displaced. Understanding such move- 
ments is important for analysing quads and for defining the 
limit curves in the next section. Fig. |5] shows the locations 
of the images, labelled by A, B, C, D, and E (middle panel), 
when the source is at the centre of the astroid caustic (left 
panel). Despite having five images, the system is called a 
quad because the central image is usually de-magnified and 
lies near the lens galaxies, making it nearly observationally 



invisible 2 . The arrival time delay contours in the right panel 
show that the image locations are at the time delay ex- 
trema or saddles, except for the extrema where the su rface 
mass densities are non-smooth jKochanek et alJ I2004T) . At 
the centroids of Gl and G2 whose locations are given in Ta- 
blcQ the time delay achieves local maxima, but there are no 
corresponding images because the surface mass densities are 
singular at the galaxy centroids in the model described by 
equation JSJ. Ignoring the central image (E, which is finitely 
de-magnified), the two images (C and D) inside the critical 
curve are time delay saddles, and the two images (A and 
B) outside the critical curve are time delay minima. This is 
true in general for quads. 

Fig. [3] shows the image locations and the time delay 
contours as the source moves across a fold from within the 
caustic. As the source approaches a fold, two of the images 
(B and C for the upper fold of interest) that are separated 
by the critical curve come together. When the source is on 
the fold, the two images merge to become one at the cor- 
responding point on the critical curve. Finally, when the 
source moves across the fold, the merged image disappears. 
The merging and disappearance of the two images can be ex- 
plained using the lemniscate time delay contour (the saddle 
with two minima) in the right panels. When the source ap- 
proaches a fold, the time delay saddle of the lemniscate joins 
with one of its two associated local minima; after the source 
crosses the fold, only one time delay minimum remains. 

Fig. |1] shows the image locations and the arrival time 
delay contours as the source moves from within the astroid 
caustic across a cusp in a direction that is roughly along the 
semi-major axis of the lens distribution. As the source ap- 
proaches the cusp, three of the images (A, B, and C in this 
case) come together. Two images (A and B) are outside and 
one image (C) is inside the critical curve. When the source 
is on the cusp, the three images become one on the criti- 
cal curve. Finally, when the source moves across the cusp, 
one image remains outside the critical curve. (We label the 
remaining image by the one that comes alphabetically first 
among the three merging images.) The time delay contours 
in the right panels depict this behaviour: the time delay sad- 
dle of a lemniscate merges simultaneously with both of its 
two minima and leaves a single minimum in the end. 

Fig.|S]is similar to Fig.0]but with the source displacing 
toward a cusp that is roughly along the semi-minor axis of 
lens distribution. The three merging images now have one 
image (B) outside and two images (C and D) inside the crit- 
ical curve (shown in middle panels). In terms of the time 
delay contours (right panels), this corresponds to the simul- 
taneous merging of the saddle of the lemniscate with one 
of its minima and with the saddle of the enclosing limagon, 
leaving only the limagon saddle in the end. 

2.2.3 Inner and outer limits 

The movements of the image l ocations shown in F igs.EltoE] 
allow us to define limit curves (Blandfo rd &i Naravanl l986"). 
Consider moving a hypothetical point source on the caustic 

2 We refer the reader to lWinn. Rusin fc W ambsganssl J2004I) for 
candidates of central image detections in gravitational lens sys- 
tems. 
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Figure 2. Lef t pane l: source is near the centre within the astroid caustic of the B1608+656 SPLEl+D(isotropic) model in 
Koopman s et alJ l200.ll) . Middle panel: the corresponding five images (A, B, C, D, and E), the lens galaxy positions (Gl and G2) 
indicated by open triangles, and the critical curves. Right panel: crucial time delay contours for demonstrating Format's principle. The 
time delay at each image position is a minimum (L for "low") or a saddle (S). The scales on the source plane and image plane are 
different due to magnification of the images. 
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Figure 3. Left panels: source position displaced across a fold from inside (top) to outside (bottom) of the astroid caustic curve of 
B1608+656 SPLEl+D(isotropic) model. Middle panels: image positions (A, B, C, D, and E) corresponding to the source positions shown 
in the left panels, lens galaxy positions (Gl and G2) indicated by open triangles, and the critical curves. Right panels: corresponding 
time delay contours. Letter L (for low) or S at each image location represents a time delay minimum or saddle, respectively. 
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Figure 4. Left panels: source position displaced across a cusp approximately along the semi-major axis from inside (top) to outside 
(bottom) of the astroid caustic curve of B1608+656 SPLEl+D(isotropic) model. Middle panels: image positions (A, B, C, D, and E) 
corresponding to the source positions shown in the left panels, lens galaxy positions (Gl and G2) indicated by open triangles, and the 
critical curves. Right panels: corresponding time delay contours. Letter L (for low) or S at each image location represents a time delay 
minimum or saddle, respectively. 



curve. As the source traces around the folds of the caustic, 
the two non-merging images trace out the limit curves. For 
the astroid, the non-merging image inside the critical curve 
is on the inner limit and the image outside the critical curve 
is on the outer limit. For the deltoids, both non-merging 
images are outside the corresponding critical curves. The 
deltoids thus have only outer limits composed of two images 
and no inner limits. Fig.|S]is the p lot of the limit cu r ves fo r 
the SPLEl-l-D(isotropic) model in iKoopmans et al.l J2003I) . 
The inner limit and outer limit for the astroid are shown 
in green and orange, respectively. The outer limit for the 
deltoids are shown in cyan. 

We focus only on the limit curves of the astroid since 
they are typical for elliptical lens mass distributions. Both 
the inner and the outer limits are tangent to the critical 
curve twice, corresponding to source placement at the cusps 
of the caustic. The limit curves mark the boundary of the 
region containing four images. 

2.2.4 Isophotal separatrices 

An isophote is an intensity contour. We assume the source 
intensity distribution has a single maximum with nested, 
non-crossing contours. An isophotal separatrix on the im- 



age plane corresponds to a source intensity contour that is 
tangent to the caustic curve. The isophotes must cross at 
the critical curve and be tangent to the limit curves as we 
explain below. 

Consider an extended elliptical source intensity distri- 
bution centred at (/? s i , /3 S 2) = (0.088,-1.069) with an axis 
ratio of 0.634 and a semi-major axis position angle of 22.1 
degrees 3 . The left panel in Fig. |7| shows four coloured inten- 
sity contours of the extended source. The two intermediate 
isophotes are very close together (light blue and dark blue). 
The right panel in Fig.|7|shows the mapped isophotes (same 
colours) with the critical curves (black) and limit curves 
(red). Each coloured set of isophotes must intersect at the 
critical curve and be tangent to the inner and outer limit. 
This is shown most clearly by the purple isophotes that con- 
sist of a lemniscate (separatrix) with two elliptical satellite 
isophotes on the image plane. The lemniscate isophote must 
cross at the critical curve, and the two satellite isophotes 
must each be tangent to either the inner or the outer limit. 



3 This source model differs from the IKoopmans et all i2003t) 
source model in the position angle, but the difference is irrele- 
vant for the purpose of describing isophotal separatrices. 
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Figure 5. Left panels: source position displaced across a cusp approximately along the semi-minor axis from inside (top) to outside 
(bottom) of the astroid caustic curve of B1608+656 SPLEl+D(isotropic) model. Middle panels: image positions (A, B, C, D, and E) 
corresponding to the source positions shown in the left panels, lens galaxy positions (Gl and G2) indicated by open triangles, and the 
critical curves. Right panels: corresponding time delay contours. Letter L (for low) or S at each image location represents a time delay 
minimum or saddle, respectively. 



To explain the crossing and tangency conditions, let us 
consider the purple isophotes in detail. The crossing point 
of the lemniscate on the critical curve corresponds to the 
tangency point of the source isophote to the astroid caustic 
curve. Recall from section 2.2.2 that two of the four images 
of a hypothetical point source merge on the critical curve as 
the source moves across the fold from within. Therefore, a 
segment of the source isophote to either side of the caustic 
tangency point will map to two segments on the image plane, 
one inside and one outside the critical curve, that connect at 
the critical curve. The entire source isophote that is within 
the caustic will thus correspond to a lemniscate crossing the 
critical curve on the image plane with one lobe inside and 
one lobe outside the critical curve. The tangencies of the im- 
age isophotes to the limit curves can be understood based 
on the definition of limit curves, which are the inner and 
outer boundaries of the four-image region that are marked 
by the two non-merging images as a hypothetical source 
traces around the caustic. The two satellite isophotes cor- 
respond to image isophotes traced by the two non-merging 
images that must touch the inner and outer limits when the 
source isophote is tangent to the caustic. Since the inner and 
outer limits are the four-image boundaries, the touchings of 
the satellite isophotes to the limit curves become tangencies. 



Similar reasoning applies to the crossings and tangencies of 
the other three sets of isophotes. 

The crossing of the isophotes at the critical curves and 
the tangency of the isophotes to the limit curves provide 
qualitative tests on how good a lens model is. 

2.2.5 Observational Data 

We use the result of section 2. 2.4 to qualitative l y test the 
SPLEl+D(isotropic) model in iKoopmans et alJ fcOO&t by 
superimposing the critical and limit curves of the model on 
the intensity contours of the observational data. Fig.|H|shows 
the isophotal separatrices (in black in various line styles) of 
the deconvolved residua l HST /F814W image of B1608+656 
jKoopmans et alJ 120031) with the critical curves (red) and 
limit curves (green, orange, cyan). We check the crossing and 
tangency conditions for each of the four sets of isophotal sep- 
aratrices, using Fig.|7|as a guide for the approximate cross- 
ing and tangency locations. For the dashed isophotes, the 
conditions for the crossing of the separatrix at the critical 
curve and the tangency to the limit curves are violated. For 
the solid isophotes, the crossing at (81,62) ~ (—0.8,-1.1) 
is not at the critical curve, but the tangency requirements 
at ~ (0.9,-1.4) and ~ (0.4,0.3) are satisfied within the 
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Figure 7. Left panel: four isophotes (in colours) of an extended source intensity distribution that are tangent to the astroid caustic 
curve (black). Right panel: the mapping of the isophotes in the left panel with the critical curve (black) and limit curves (red). These 
isophotes must cross at the critical curve and their satellite isophotes must be tangent to the limit curves. 
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Figure 6. The limit curves are plotted with the critical curves 
(in black) for the SPLEl+D(isotropic) model of B1608+656. The 
orange (green) curve is the outer (inner) limit associated with the 
astroid. The limit curves are each tangent to the critical curve of 
the astroid twice. The cyan curves are the outer limits of the 
deltoids. 

noise. For the dotted isophotes, the crossing at ~ (0.7, —1.9) 
is at the critical curve within the noise, but the isophotes 
near ~ (—0.5,-0.9) and ~ (1.1,0.2) are not tangent to 
the limit curves. Lastly, for the long-dashed isophotes, the 
crossing at ~ (1.3, —0.6) is on the critical curve, and the 



isophotes near ~ (—0.5, —0.6) and ~ (—0.3, —2.4) are tan- 
gent to the limit curves, within the nois e. Therefore, the 
SPLE 1+D (isotropic) model proposed by iKoopmans et all 
(2003) satisfies the crossing and tangency conditions stated 
in section 2.2.4 for some, but not all, of the isophotal separa- 
trices. As a result, the model proposed bv lKoopmans et alJ 
( 2003) must not represent the true lens potential of the sys- 
tem, especially in the regions where the crossings and tan- 
gencies fail. Recall that we need an accurate model of the 
lens potential to calculate the Hubble constant. In the next 
section, we examine a method of potential correction. 



3 POTENTIAL RECONSTRUCTION 

3.1 Theory of potential reconstruction 

Th e method of potential r econstruction was first suggested 
bv iBlandford et al.l (1200 ll) . Following the notation in sec- 
tion 2.1, let 1(9) be the observed image intensity of a grav- 
itational lens system with an extended source. For a given 
potential model, tp(0) , one can obtain the b est-fitting source 
intensity distribution dWarren fc Dvei2003t) . Let 7Q3) be the 
source intensity translated to the image plane via the poten- 
tial model, ip(0). We define the intensity deficit on the image 
plane by 

51(6) =1(9) -I ((3). (10) 

The intensity deficit is zero everywhere with the true lens 
potential distribution. 

Consider a lens potential model that is perturbed from 
the true potential, ipo(6), by S-if)(d): 

iP(6) = M0) + H(0)- (11) 



8 S. H. Suyu and R. D. Blandford 




-3 

-2-10 1 2 

01 (arcsec) 

Figure 8. The deconvolved resi dual ffi>T/F814W image of 
B1608+656 (Koopmans et 21 II I.J'). The isophotal separatrices 
(in black in various line styles) are shown with the critical 
curves (red) and limit curves (green, orange, and cya n) of the 
SPLEl+D(isotropic) model in lKoopmans et alj J2003I) . Some of 
the isophotal separatrices are not intersecting at the critical curve 
of the model and some of the satellite isophotes are not tangent 
to the limit curves of the model. 



We can correct the potential model perturbatively by solving 
for the perturbation 5ip(&). For a given image (fixed and 
1(0)), we can relate a change in position on the source plane, 
5(3, to the potential perturbation using the lens equation Q: 



50 = - 



85jj(0) 
89 ' 



(12) 



Expanding 1(0) to first order in 5/3 and using equation l|12l) 
in equation I10L we obtain 



51(0) 



dm 

d/3 



■ 5(3 



81(0) 854,(0) 



0/3 



80 



(13) 



The source intensity gradient dI g^ implicitly depends on 
the potential model ip(0) since the source position (3 (where 
the gradient is evaluated) is related to ip(0) via the lens 
equation To first order, using the perturbed model ip(0) 
is equivalent to using the true model ipo(0) in the evaluation 
of the source intensity gradient dI g^ ■ 

We can solve equation 1131 for the potential correction, 
5tp(0), provided that we start at a potential model that is 
close to the true potential. (We quantify what "close" means 
in the next section.) One method to solve for the potential 
correction is to integrate along the characteristics of the par- 
tial differential equation 11311 . The solution is 



54,(0) =5^(0 A ) + 



d6 s 5I(0) 

I dI(/3) p 

°a | a/3 | 



where 



d6 s = [del + de l) 



2\l/2 



(14) 



(15) 



dI(/3) 



8/3 



(8I(0)\ (8I((3)\ 
\ d(3 x ) + \ 8(3 2 J 



(16) 



and a is an arbitrary reference point that is conveniently 
chosen to be at the location of one of the images, say A. 
(The reference point is arbitrary because the potential is 
determined up to a constant.) The characteristic curves, on 
which we must integrate to obtain the potential correction, 
are given by curves that satisfy 

d6i _ 81/8(3! 

~dT 2 ~ 



(17) 



81/8(32 

Each point on a characteristic curve thus follows the source 
intensity gradient (evaluated at the corresponding source lo- 
cation given by the lens equation Q) that is directly trans- 
lated to the image plane without distortions via the magni- 
fication matrix. Due to the direct translation of the source 
intensity gradient, the characteristic curves differ from the 
curves on the image plane that map to the source intensity 
gradient curves. The structure of the characteristic curves 
allows us to determine whether the potential solution given 
by equation 1141 is unique. This is demonstrated in the next 
section. 

We can repeat the process for a perturbative and itera- 
tive potential reconstruction method. We expect the poten- 
tial to be closer to the true potential after each iteration, 
which is indicated by a decrease in the magnitude of the 
intensity deficit. 

The potential reconstruction method is non-parametric. 
We can pixellate the potential distribution to match the 
observed image pixellation, and the potential correction at 
each pixel is given by equation 1141 1. 

To summarise, the four steps for the method are: (i) 
start with a potential model close to the true potential, (ii) 
calculate the intensity deficit (equation 1101 ) of each pixel, 
(iii) calculate the potential correction of each pixel (equa- 
tion 1141 ) by integrating along the characteristics (equation 
117> ). (iv) obtain the corrected potential and repeat the pro- 
cess (steps (ii) to (iv)) until the intensity deficit approaches 
zero. In the next section, we examine a quadruply imaged 
toy model to test the method of potential reconstruction. 

3.2 Example Toy Model 

To demonstrate the method of potential reconstruction dis- 
cussed in the previous section, we consider a toy model 
with a simple lens potential that produces a quad like 
B1608+656. 

The toy system has a non-singular isothermal ellipsoid 
lens whose potential takes the form: 

1/2 



4,(e 1 ,e 2 ) = (0? + 20I + O.1 



(18) 



We take the perturbed potential to be the original potential 
that is rotated clockwise by 1.1 degrees. The source intensity 
distribution has elliptical contours with axis ratio of 0.634 
and position angle of 147.2 degrees. The source nucleus is 
located at ((3is,02s) = (0.1,0.05) and has an intensity peak 
of 100, in arbitrary units. We assume the data is perfect 
with no noise, but we discretize the image plane region [- 
2,2] x [-2,2] into a 201x201 grid in order to correct for the 
perturbation of every pixel. In Fig. [5] the left panel shows 
the caustic curves (dashed) of the original potential and the 
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Figure 10. The time delay contour associated with the source 
nuclear position of the toy model. The image locations of the 
source nucleus (see Fig. fright panel) are at time delay saddles 
(S), minima (L) or maxima (H). 

source intensity contours (dotted), and the right panel shows 
the corresponding critical curves (dashed) and image inten- 
sity contours (dotted). Analogous to B1608+656, there is 
an astroid caustic in the left panel. The additional elliptical 
caustic curve is due to the non-singular nature of the lens 
potential. Different regions separated by the caustic curves 
have different image multiplicities. In the enclosed region 
intersected by the astroid and elliptical caustic curves, a 
source has five images on the image plane. In the region 
within the caustic curves excluding the intersection, a source 
has three images. In the region outside the caustic curves, 
a source has one image. The astroid caustic is mapped to 
the outer critical curve and the elliptical caustic is mapped 
to the inner critical curve. As for BI608+656, we focus on 
the astroid caustic and the outer critical curve. Among the 
isophotes in the right panel, the four isophotal separatrices 
that are shown match to the four isophotes tangent to the 
astroid caustic in the left panel. The separatrices intersect at 
the outer critical curve, as required (section 2.2.4). Fig. 1101 
shows the arrival time delay contour of the source nucleus 
of the toy model. The quad has similar time delay extrema 
(two saddles within the critical curve and two minima out- 
side the critical curve) to the SPLEl+D(isotropic) model of 
B1608+656. 

We simplify the potential correction method by using 
the original source intensity distribution and the character- 
istic fields of the original potential (instead of reconstructing 
from the perturbed potential) . In reality, we would have to 
use the reconstructed source dWarren fc Dvell2003D and the 
characteristic fields of the perturbed potential. This would 
involve simultaneously determining the source and lens po- 
tential distributions and investigating the partial degeneracy 
between them, which are beyond the scope of this paper. We 
use the simplifying assumptions on the source intensity and 




Figure 11. The characteristic fields of the toy potential model. 
The attractors are associated with images that are time delay 
minima/maxima and the repellors are associated with time delay 
saddles. 



characteristic curves as the first step to testing the method 
of potential reconstruction via integration along character- 
istics. Only if the method works robustly in this simplified 
regime is the consideration of the more general problem rel- 
evant. 

Fig. II II shows the characteristic field given by equation 
1171 1. The field has "attractors" (where field lines come to- 
gether) and "repellors" (where field lines curve away) at the 
image locations of the source nucleus. Using equation J1J and 
noting that the Jacobian matrix of T(0,/3) with respect to 
is equivalent to A in equation @ up to a constant coef- 
ficient, one can show that the attractors (or repellors) are 
associated to time delay minima/maxima (or saddles) for a 
source distribution that has non-crossing intensity contours. 
A comparison between Fig. llUl and Fig. llll confirms this fact. 

We need to follow along the characteristics to correct 
for the potential perturbation given by equation 111411 . In 
Fig. 1111 almost all of the characteristic curves end at one 
of the three attractors; but there are special characteristic 
curves that connect the attractors and repellors. These four 
connecting characteristics between the four images (exclud- 
ing the central image), shown in the right panel of Fig. [^] in 
solid lines, allow us to fix the potential offsets between the 
images and hence uniquely determine the potential up to a 
constant. The left panel of Fig.[5]shows the mapping of these 
connecting characteristics onto the source plane (solid lines) . 
As one may expect, the mapping of each of the connecting 
characteristics between an attractor and a repellor is a loop 
on the source plane that is tangent to the astroid caustic 
curve due to the connecting characteristics intersecting the 
outer critical curve. 

In addition to the characteristic curves, the intensity 
deficit is required for potential correction in equation 11411 . 
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Figure 9. Left: the caustic curve (dashed) of the original toy potential model with the intensity contours (dotted) of the source. Four 
of the intensity contours are tangent to the caustic curves. The four mappings of the connecting characteristics (solid) are each tangent 
to the caustics. Right: the critical curves (dashed) of the original toy potential model and the image intensity contours (dotted), four 
of which are isophotal separatrices intersecting at the outer critical curve. The four connecting characteristics (solid) between the four 
images each cross the outer critical curve once. 



To get the intensity deficit defined in equation 1101 for the 
pixels on the image plane, first we use the perturbed po- 
tential model, the lens equation £Q, an< i the original source 
intensity distribution to get 7(/3), then we subtract it from 
1(9) obtained from the original potential. Fig. H2l shows the 
initial intensity deficit and the initial potential perturbation 
(8ip(6) in equation mil l before the perturbative and iter- 
ative potential correction, in the top left and bottom left 
panels, respectively. We use plots of 5ip(9) to check that the 
perturbation approaches zero after corrections. 

In each potential reconstruction iteration, we use the 
current perturbed potential model to obtain the intensity 
deficit (51(0)) and the source intensity gradient ( | ^jjp | ) a ^ 
every pixel on the image plane; we then use equation 1141 
to correct the perturbed potential by integrating along the 
characteristic curves of the original potential model. Two 
iterations are performed and the resulting intensity deficit 
and potential perturbation after each iteration are shown 
in Fig. 1121 The middle and right panels show the intensity 
deficit (potential perturbation) in the top (bottom) after 1 
and 2 iterations, respectively. The middle and right panels 
are plotted on the same scales as that in the left panels. 
Comparing the right panels to the left panels, the inten- 
sity deficit and potential perturbation converge to zero after 
two iterations (apart from numerical error), signifying that 
the method of potential reconstruction along characteristics 
works in theory with perfect data. 

A possible limitation to this method is that the inten- 
sity deficit needs to be zero at the image locations; other- 
wise, according to equation 1141 . the integrand diverges at 
the image locations, which are the end points of integration. 
For the example above, we are saved from this divergence 
by discretizing the image plane and thus only reaching the 



image points within some tolerance, but never ending at the 
image (divergent) points. The potential correction is most 
significant near the image points for any non-zero intensity 
deficit in the region. Therefore, integrating along the charac- 
teristics may place limitations on the magnitude of potential 
perturbation that we can correct, which we discuss in the 
next section. 

This method of potential reconstruction works only for 
small potential perturbations like the example we consid- 
ered where the perturbation magnitude is on average (over 
the image grid) 0.5 per cent of the original potential. By 
increasing the rotation of the original potential distribution 
to get the perturbed potential (that is, increasing the per- 
turbation), we require more iterations for convergence, as 
expected. When the rotation the original potential gets to 
~ 4.5 degrees, which corresponds to an average potential 
perturbation magnitude of ~ 1.5 per cent, the method ceases 
to converge. Therefore, the method of potential correction 
by integrating equation 1141 along characteristics works in 
theory with perfect data with a small (< 1 per cent) po- 
tential model error. Unless a better algorithm is found for 
treating larger potential pertu rbation and real data with 
noise, the method proposed bv lBlandford et al.l (l200lT) will 
not in practice be useful. 

The example toy model considered provides a practical 
insight into the theory of potential reconstruction. In reality, 
we do not have useful data everywhere due to the presence of 
noise; for an extended source, we can observe emission in an 
Einstein ring connecting the four images. Based on the anal- 
ysis of this section, the Einstein ring must be large enough 
to enclose the connecting characteristics in order to obtain 
proper potential offsets between the images. This condition 
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Figure 12. Top row, left to right: the intensity deficit before potential correction, after 1 iteration and after 2 iterations of potential 
correction. The maximum initial intensity deficit (top left) is 14 near the image positions (the peak nuclear source intensity is 100). 
Bottom row, left to right: potential perturbation before correction, after 1 iteration and after 2 iterations of correction. The initial potential 
perturbation magnitude (bottom left) is on average around 0.5 per cent of the original potential. Since the potential is determined up to 
an arbitrary constant, the potential perturbation is plotted with respect to the mean to enhance small scale features. The plotting scales 
of the middle and right panels (after corrections) are the same as the left panels (before corrections) for comparison. 



must hold for any potential reconstruction algorithm based 
on equation (15). 



4 CONCLUSIONS AND FURTHER WORK 

We have considered the gravitational lens system 
B1608+656 to investigate the properties of quads. We 
have defined limit curves as the loci of non-merging images 
as the source traces the caustic curve. For the typical 
astroid caustic curve of a quad, there are inner and outer 
limits (relative to the critical curve) that are each tangent 
to the critical curve twice. We have shown that isophotes 
that are tangent to the astroid caustic curve on the source 
plane map to isophotal separatrices on the image plane. 
These separatrices must intersect on the critical curve 
and their associated satellite isophotes must be tangent to 
the limit cu rves. We have sh o wn th at the current model 
proposed by iKoopmans et alJ fcOO&t for B1608+656 does 
not satisfy these qualitative constraints for some of the 
isophotal separatrices. 

We have investigated a perturbative and itera- 
tive method of pot ential reconstruction proposed by 
iBlandford et alJ (1200 ll) . The method requires solving a par- 
tial differential equation for the potential correction, which 
we have done by integrating along the characteristic curves. 



We have used a toy model that is a quad like B1608+656 
to test the method, assuming perfect data. For small per- 
turbations whose magnitudes are on average < 1 per cent of 
the original potential, the method has worked and we have 
had the perturbed potential converging to the true poten- 
tial. However, the method has failed to converge when the 
perturbation magnitude increases to around 1.5 per cent of 
the original potential. This may be due to the non-zero in- 
tensity deficit near the image locations which restricts the 
integration along characteristics. We hope to use the knowl- 
edge we have acquired about the anatomy of the quads and 
the characteristic fields of the potential correction equation 
to find a more robust method of potential correction that 
can be applied to real data with noise. 
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